# This script performs link-predictive analysis of the PE network
# Last edited: 9/24/2013
# Edited by: Bruce Desmarais

library(oc)

# Read in network data
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/netList.RData")
# amat is the adjacency matrix recording the number of events both attend
# vdat is the vertex dataset containing the variables in DSMKdata

### Cosponsorship ###
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/netListC.RData")

load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/RollCallRes.RData")
estsL <- ests
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/RollCallResNoLog.RData")
estsNoLog <- ests
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/RollCallResNL.RData")
estsNL <- ests
load("~/Dropbox/professional/Research/Active/SenPE/Data/cleanData/RollCallResNLNoLog.RData")
estsNLNoLog <- ests

# phi coefficient

congs <- c(96:105)

for(i in congs){
filei <- paste("~/Dropbox/professional/Research/Active/SenPE/Tex/RCpredict_",i,".pdf",sep="")
fileNLi <- paste("~/Dropbox/professional/Research/Active/SenPE/Tex/RCpredictNL_",i,".pdf",sep="")

penet <- netList[[i-95]]$amat

xpe <- penet[lower.tri(penet)]
xsp <- seq(min(xpe),max(xpe),length=1000)

betaL <- coef(estsL[[i-95]])[2] 
betaNoLog <- coef(estsNoLog[[i-95]])[2] 
betaNL <- coef(estsNL[[i-95]])[2] 
betaNLNoLog <- coef(estsNLNoLog[[i-95]])[2] 

prL <- betaL*log(xsp+1) - mean(betaL*log(xsp+1)) + mean(estsL[[i-95]]$fitted.values)
prNoLog <- betaNoLog*xsp  - mean(betaNoLog*xsp)+mean(estsNoLog[[i-95]]$fitted.values)
prNL <- plogis(betaNL*log(xsp+1) - mean(betaNL*log(xsp+1)) + mean(estsNL[[i-95]]$fitted.values))*2-1
prNLNoLog <- plogis(betaNLNoLog*xsp - mean(betaNLNoLog*xsp) + mean(estsNL[[i-95]]$fitted.values))*2-1


sigL <- summary(estsL[[i-95]])$pgreq[2] < .05
sigNoLog <- summary(estsNoLog[[i-95]])$pgreq[2] < .05
sigNL <- summary(estsNL[[i-95]])$pgreq[2] < .05
sigNLNoLog <- summary(estsNLNoLog[[i-95]])$pgreq[2] < .05

yl <- range(c(prL,prNoLog,prNL,prNLNoLog))

pdf(filei,height=1.6, width=5.75, family="Times", pointsize=12.25)
par(las=1,mar=c(2.5,5.5,.75,.5),cex.lab=1.5,cex.axis=1.15)
plot(xsp,prL,lty=1,col=ifelse(sigL,"black","grey50"),lwd=4,type="l",xlab="",ylab="",ylim=yl)
lines(xsp,prNoLog,lty=2,col=ifelse(sigNoLog,"black","grey50"),lwd=4)
#rug(xpe)
#title(xlab="Joint Press Events",line=2.25)
title(ylab="Corr",line=4)
dev.off()

pdf(fileNLi,height=1.6, width=5.75, family="Times", pointsize=12.25)
par(las=1,mar=c(2.5,5.5,.75,.5),cex.lab=1.5,cex.axis=1.15)
plot(xsp,prNL,lty=1,col=ifelse(sigNL,"black","grey50"),lwd=4,type="l",xlab="",ylab="",ylim=yl)
lines(xsp,prNLNoLog,lty=2,col=ifelse(sigNLNoLog,"black","grey50"),lwd=4)
#rug(xpe)
#title(xlab="Joint Press Events",line=2.25)
title(ylab="Corr",line=4)
dev.off()


}

library(xtable)

effTab <- NULL
for(i in 1:length(estsL)){
	coefsi <- round(coef(estsL[[i]])[3:7],dig=4) 
	signifi <- estsL[[i]]$pgreqabs[3:7]
	coefsi <- paste(coefsi,ifelse(signifi < .1,"*",""),sep="")
	effTab <- rbind(effTab,coefsi)
	
}




